Il dataset analizzato è “Descrizione tracciato record dati comunali” dell’ISTAT sul numero di decessi in Italia per i mesi da Gennaio ad Aprile e per gli anni dal 2015 al 2020.
La prima modifica effettuata al dataset è stato unpivot dei dati per la classe sesso_anno. Sucessivamente le informazioni riguardo al genere e all’anno di raccolta del dato sono state separate in due differenti colonne, sono state create le colonne DATA, MESE e GIORNO. Infine sono state eliminate le colonne non più necessarie.
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyr)
library(ggplot2)
library(lubridate)
##
## Attaching package: 'lubridate'
## The following objects are masked from 'package:dplyr':
##
## intersect, setdiff, union
## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, union
library(forecast)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
library(zoo)
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(vctrs)
library(fitdistrplus)
## Loading required package: MASS
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
## Loading required package: survival
## Loading required package: npsurv
## Loading required package: lsei
library(actuar)
##
## Attaching package: 'actuar'
## The following object is masked from 'package:grDevices':
##
## cm
df_path <- '/Volumes/HDD_Ale/DSLab/dati-giornalieri-comune/comune_giorno.csv'
data_path_pop_tot <- "/Volumes/HDD_Ale/ds_lab/dataset/tot_pop.csv"
# df_path <- 'C:\\Users\\stefa\\Desktop\\dsLab\\local\\dataset\\comune_giorno.csv'
# data_path_pop_tot <- "C:\\Users\\stefa\\Desktop\\dsLab\\local\\dataset\\tot_pop.csv"
df <- read.csv(df_path)
df_pop <- read.csv(data_path_pop_tot)
# sistemo il dataset in forma più comoda
df <- df %>%
gather(key = "SESSO_ANNO", value = "DECESSI", MASCHI_15:FEMMINE_20)
df %>% separate(SESSO_ANNO, c("SESSO", "ANNO"), "_") -> df
df %>% mutate(DATA = as.Date(paste0("0", GE, "20", ANNO), format = "%m%d%Y")) -> df
df$GIORNO <- as.numeric(substr(as.character(df$GE), 2, 3))
df$MESE <- as.numeric(substr(as.character(df$GE), 1, 1))
df$GE <- NULL
df$TOTALE_15 <- NULL
df$TOTALE_16 <- NULL
df$TOTALE_17 <- NULL
df$TOTALE_18 <- NULL
df$TOTALE_19 <- NULL
df$TOTALE_20 <- NULL
summary(df)
L’analisi dei decessi è focalizzata sulla regione Lombardia.
La prima analisi effettuata per dimostrare che c’è stato un incremento significativo dei decessi nel 2020 rispetto agli anni precedenti è un’analisi della varianza tra i gruppi anno 15, 16, 17, 18, 19, 20.
I dati relativi ad Aprile sono stati inoltre eliminati in quanto non totalmente disponibili nell’anno 2020.
La colonna DECESSI riporta diversi valori pari a 9999. Questo dato è un valore dummy per indicare la mancanza del dato in fase di raccolta. Essendo che vogliamo dimostrare la differenza in media delle distribuzioni tra i gruppi consideriamo solamente i comuni che non presentano questa anomalia.
Il dataset così pulito è formato da 622 comuni Lombardi.
df %>% filter(NOME_REGIONE == 'Lombardia',
DECESSI != 0,
MESE %in% c(1,2,3))%>%
dplyr::select(NOME_COMUNE,ANNO,MESE,DECESSI,DATA) -> df_comuni
#summary(df_comuni)
df_comuni[df_comuni$DECESSI == 9999,] %>% dplyr::select(NOME_COMUNE) %>% distinct -> comuni_9999
df_comuni_NO_9999 <- subset(df_comuni, !(df_comuni$NOME_COMUNE %in% comuni_9999$NOME_COMUNE))
summary(df_comuni_NO_9999)
## NOME_COMUNE ANNO MESE DECESSI
## Milano : 8522 Length:110978 Min. :1.000 Min. : 1.000
## Brescia: 2996 Class :character 1st Qu.:1.000 1st Qu.: 1.000
## Bergamo: 2196 Mode :character Median :2.000 Median : 1.000
## Monza : 1916 Mean :2.038 Mean : 1.188
## Como : 1511 3rd Qu.:3.000 3rd Qu.: 1.000
## Cremona: 1467 Max. :3.000 Max. :17.000
## (Other):92370
## DATA
## Min. :2015-01-01
## 1st Qu.:2016-02-29
## Median :2018-01-24
## Mean :2017-11-03
## 3rd Qu.:2019-03-21
## Max. :2020-03-31
##
#622 comuni senza 9999
length(unique(df_comuni_NO_9999$NOME_COMUNE))
## [1] 622
Aggregando per comune ed anno, sommando il numero di decessi scopriamo la presenza di forti outliers. Questo suggerisce di differenziare l’analisi in categorie differenti, in modo da evitare che le code pesino eccessivamente sulle distribuzioni.
#Tutti mesi (1,2,3)
#61.5 = 30 + (30-9)*1.5
df_aggr <- plyr::ddply(df_comuni_NO_9999, c("NOME_COMUNE","ANNO"), function(x) colSums(x[c("DECESSI")]))
summary(df_aggr)
## NOME_COMUNE ANNO DECESSI
## Acquafredda : 6 Length:3726 Min. : 1.00
## Acquanegra Cremonese : 6 Class :character 1st Qu.: 9.00
## Acquanegra sul Chiese: 6 Mode :character Median : 15.00
## Adrara San Martino : 6 Mean : 35.39
## Agnadello : 6 3rd Qu.: 30.00
## Agrate Brianza : 6 Max. :4224.00
## (Other) :3690
Per far ciò si considera quindi la classe dei comuni con bassa mortalità, utilizzando come punto di riferimento il massimo \(Q_3 + 1.5\cdot(Q_3-Q_1) = 61.5\).
df_aggr[df_aggr$DECESSI > 61,] %>% dplyr::select(NOME_COMUNE) %>% distinct -> comuni_alti
df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_alti$NOME_COMUNE),]$NOME_COMUNE -> comuni_bassi
df_aggr_bassi <- subset(df_aggr, df_aggr$NOME_COMUNE %in% comuni_bassi)
summary(df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_bassi),])
#188.4 = 98.75 + (98.75-39)*1.5
df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_bassi) & df_aggr$DECESSI > 188,] %>% dplyr::select(NOME_COMUNE) %>% distinct -> comuni_alti
df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_alti$NOME_COMUNE) & !(df_aggr$NOME_COMUNE %in% comuni_bassi),]$NOME_COMUNE -> comuni_medi
df_aggr_medi <- subset(df_aggr, df_aggr$NOME_COMUNE %in% comuni_medi)
df_aggr_alti <- subset(df_aggr, df_aggr$NOME_COMUNE %in% comuni_alti$NOME_COMUNE)
Essendo ancora elevata la presenza di outliers si procede con la classe dei comuni con media mortalità, utilizzando come punto di riferimento il massimo \(Q_3 + 1.5(Q_3-Q_1) = 188.4\). I restanti comuni formeranno la classe ad alta mortalità. Infine i comuni di Milano, Bergamo e Brescia verranno considerati a parte in quanto il numero di decessi è molto più elevato della distribuzione dell’ultima classe costruita.
par(mfrow=c(1,1))
boxplot(DECESSI~ANNO,data=df_aggr_bassi, outline = TRUE, main="# Decessi nei Comuni Lombardi con Basso numero di morti",
xlab="Anni", ylab="# Decessi", col="lightblue")
boxplot(DECESSI~ANNO,data=df_aggr_medi, outline = TRUE, main="# Decessi nei Comuni Lombardi con Medio numero di Morti",
xlab="Anni", ylab="# Decessi", col="lightblue")
boxplot(DECESSI~ANNO,data=df_aggr_alti[!(df_aggr_alti$NOME_COMUNE %in% c("Milano","Bergamo","Brescia")),], outline = TRUE, main="# Decessi nei Comuni Lombardi con Alto numero di Morti",
xlab="Anni", ylab="# Decessi", col="lightblue")
boxplot(DECESSI~ANNO,data=df_aggr_alti[df_aggr_alti$NOME_COMUNE %in% c("Milano","Bergamo","Brescia"),], outline = TRUE, main="# Decessi nel Comune di Milano, Bergamo e Brescia",
xlab="Anni", ylab="# Decessi", col="lightblue")
Dall’analisi grafica sui boxplot possiamo notare la presenza di un’incremento nell’anno 2020. Per approfondire le cause dell’aumento del numero di decessi si osservano i dati mensili separatamente, applicando la logica delle classi di mortalità.
#Gennaio----------------------------------------------------------------------------------------------------------------------------
df_aggr <- plyr::ddply(df_comuni_NO_9999[df_comuni_NO_9999$MESE == 1,], c("NOME_COMUNE","ANNO"), function(x) colSums(x[c("DECESSI")]))
#summary(df_aggr)
#20.5 = 10 + (10-3)*1.5
df_aggr[df_aggr$DECESSI > 20,] %>% dplyr::select(NOME_COMUNE) %>% distinct -> comuni_alti
df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_alti$NOME_COMUNE),]$NOME_COMUNE -> comuni_bassi
df_aggr_bassi <- subset(df_aggr, df_aggr$NOME_COMUNE %in% comuni_bassi)
#summary(df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_bassi),])
#70.5 = 39 + (39-18)*1.5
df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_bassi) & df_aggr$DECESSI > 70,] %>% dplyr::select(NOME_COMUNE) %>% distinct -> comuni_alti
df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_alti$NOME_COMUNE) & !(df_aggr$NOME_COMUNE %in% comuni_bassi),]$NOME_COMUNE -> comuni_medi
df_aggr_medi <- subset(df_aggr, df_aggr$NOME_COMUNE %in% comuni_medi)
df_aggr_alti <- subset(df_aggr, df_aggr$NOME_COMUNE %in% comuni_alti$NOME_COMUNE)
par(mfrow=c(1,1))
boxplot(DECESSI~ANNO,data=df_aggr_bassi, outline = TRUE, main="# Decessi nei Comuni Lombardi con Basso numero di morti",
xlab="Anni", ylab="# Decessi", col="lightblue")
boxplot(DECESSI~ANNO,data=df_aggr_medi, outline = TRUE, main="# Decessi nei Comuni Lombardi con Medio numero di Morti",
xlab="Anni", ylab="# Decessi", col="lightblue")
boxplot(DECESSI~ANNO,data=df_aggr_alti[df_aggr_alti$NOME_COMUNE != "Milano",], outline = TRUE, main="# Decessi nei Comuni Lombardi con Alto numero di Morti",
xlab="Anni", ylab="# Decessi", col="lightblue")
boxplot(DECESSI~ANNO,data=df_aggr_alti[df_aggr_alti$NOME_COMUNE == "Milano",], outline = TRUE, main="# Decessi nel Comune di Milano",
xlab="Anni", ylab="# Decessi", col="lightblue")
Possiamo notare che graficamente non risultano differenze tra le distribuzioni. In particolare risulta un incremento, di cui andrebbe indagata la significatività nel 2017.
#Febbraio----------------------------------------------------------------------------------------------------------------------------
df_aggr <- plyr::ddply(df_comuni_NO_9999[df_comuni_NO_9999$MESE == 2,], c("NOME_COMUNE","ANNO"), function(x) colSums(x[c("DECESSI")]))
#summary(df_aggr)
#18 = 9 + (9-3)*1.5
df_aggr[df_aggr$DECESSI > 18,] %>% dplyr::select(NOME_COMUNE) %>% distinct -> comuni_alti
df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_alti$NOME_COMUNE),]$NOME_COMUNE -> comuni_bassi
df_aggr_bassi <- subset(df_aggr, df_aggr$NOME_COMUNE %in% comuni_bassi)
#summary(df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_bassi),])
#60 = 33 + (33-15)*1.5
df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_bassi) & df_aggr$DECESSI > 60,] %>% dplyr::select(NOME_COMUNE) %>% distinct -> comuni_alti
df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_alti$NOME_COMUNE) & !(df_aggr$NOME_COMUNE %in% comuni_bassi),]$NOME_COMUNE -> comuni_medi
df_aggr_medi <- subset(df_aggr, df_aggr$NOME_COMUNE %in% comuni_medi)
df_aggr_alti <- subset(df_aggr, df_aggr$NOME_COMUNE %in% comuni_alti$NOME_COMUNE)
par(mfrow=c(2,2))
boxplot(DECESSI~ANNO,data=df_aggr_bassi, outline = TRUE, main="# Decessi nei Comuni Lombardi con Basso numero di morti",
xlab="Anni", ylab="# Decessi", col="lightblue")
boxplot(DECESSI~ANNO,data=df_aggr_medi, outline = TRUE, main="# Decessi nei Comuni Lombardi con Medio numero di Morti",
xlab="Anni", ylab="# Decessi", col="lightblue")
boxplot(DECESSI~ANNO,data=df_aggr_alti[df_aggr_alti$NOME_COMUNE != "Milano",], outline = TRUE, main="# Decessi nei Comuni Lombardi con Alto numero di Morti",
xlab="Anni", ylab="# Decessi", col="lightblue")
boxplot(DECESSI~ANNO,data=df_aggr_alti[df_aggr_alti$NOME_COMUNE == "Milano",], outline = TRUE, main="# Decessi nel Comune di Milano",
xlab="Anni", ylab="# Decessi", col="lightblue")
Anche per Febbraio non vediamo particolari differenze tra il 2020 e gli altri anni. Ci aspettiamo dunque che tutta la differenza notata nelle distribuzioni totali sia apportata dal mese di Marzo.
#Marzo----------------------------------------------------------------------------------------------------------------------------
df_aggr <- plyr::ddply(df_comuni_NO_9999[df_comuni_NO_9999$MESE == 3,], c("NOME_COMUNE","ANNO"), function(x) colSums(x[c("DECESSI")]))
#summary(df_aggr)
#25,5 = 12 + (12-3)*1.5
df_aggr[df_aggr$DECESSI > 25,] %>% dplyr::select(NOME_COMUNE) %>% distinct -> comuni_alti
df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_alti$NOME_COMUNE),]$NOME_COMUNE -> comuni_bassi
df_aggr_bassi <- subset(df_aggr, df_aggr$NOME_COMUNE %in% comuni_bassi)
#summary(df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_bassi),])
#69,5 = 32 + (32-7)*1.5
df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_bassi) & df_aggr$DECESSI > 69,] %>% dplyr::select(NOME_COMUNE) %>% distinct -> comuni_alti
df_aggr[!(df_aggr$NOME_COMUNE %in% comuni_alti$NOME_COMUNE) & !(df_aggr$NOME_COMUNE %in% comuni_bassi),]$NOME_COMUNE -> comuni_medi
df_aggr_medi <- subset(df_aggr, df_aggr$NOME_COMUNE %in% comuni_medi)
df_aggr_alti <- subset(df_aggr, df_aggr$NOME_COMUNE %in% comuni_alti$NOME_COMUNE)
par(mfrow=c(2,2))
boxplot(DECESSI~ANNO,data=df_aggr_bassi, outline = TRUE, main="# Decessi nei Comuni Lombardi con Basso numero di morti",
xlab="Anni", ylab="# Decessi", col="lightblue")
boxplot(DECESSI~ANNO,data=df_aggr_medi, outline = TRUE, main="# Decessi nei Comuni Lombardi con Medio numero di Morti",
xlab="Anni", ylab="# Decessi", col="lightblue")
boxplot(DECESSI~ANNO,data=df_aggr_alti[!(df_aggr_alti$NOME_COMUNE %in% c("Milano","Bergamo","Brescia")),], outline = TRUE, main="# Decessi nei Comuni Lombardi con Alto numero di Morti",
xlab="Anni", ylab="# Decessi", col="lightblue")
boxplot(DECESSI~ANNO,data=df_aggr_alti[df_aggr_alti$NOME_COMUNE %in% c("Milano","Bergamo","Brescia"),], outline = TRUE, main="# Decessi nel Comune di Milano, Bergamo e Brescia",
xlab="Anni", ylab="# Decessi", col="lightblue")
Come intuito precedentemente il mese di Marzo è quello che più influisce sulla differenza tra gli anni nel numero di decessi. Quindi l’analisi si focalizzerà sui dati relativi a Marzo.
Ricordando le ipotesi per il test dell’ANOVA: \[
\begin{itemize}
\item[1] I soggetti sono tra loro indipendenti
\item[2] Normalità della distribuzione
\item[3] Omoschedasticità
\end{itemize}
\] La 1 è dimostrata per definizione in quanto i soggetti deceduti non verranno ri-analizzati.
L’ipotesi 2 è violata. Già dall’istogramma si nota che le distribuzioni sono non normali.
Per identificare la distribuzione del numero di decessi nei comuni lombardi è stata eseguita una statistica Goodness of Fit con stima della massima verosomiglianza dei parametri. Sono state considerate le distribuzioni teoriche seguenti: \[ \begin{itemize} \item[1] Weibull \item[2] Log Normale \item[3] Gamma \item[1] Log Logistic \item[2] Pareto \end{itemize} \]
#hp:
# 1) soggetti sono tra loro indipendenti.
# 2) normalità della distribuzione (shapiro.test(x) #pvalue > 0.05)
# 3) omoschedasticità (bartlett.test(dati, gruppi #pvalue>0.5)
#df_aggr_bassi
# 1) Rispettata perchè la persona deceduta non viene più considerata
# 2) Le distribuzioni sono chiaramente NON normali.
gos_stat <- function(x, anno){
fw <- fitdist(x, "weibull")
fln <- fitdist(x, "lnorm")
fg <- fitdist(x, "gamma")
fll <- fitdist(x, "llogis",start = list(shape = 1, scale = 500))
fp <- fitdist(x, "pareto",start = list(shape = 1, scale = 500))
#fb <- fitdist(x, "burr",start = list(shape1 = 0.3, shape2 = 1, rate = 1))
denscomp(list(fw, fln, fg, fll, fp), main = paste0("Histogram and theoretical densities, year 20",anno, sep=""))
print(fln)
return(gofstat(list(fw, fln, fg, fll, fp), fitnames = c("Weibull", "Lognormal", "Gamma", "Loglogistic", "Pareto")))
}
par(mfrow = c(3, 2))
for(i in 15:20) {
x <- df_aggr_bassi[df_aggr_bassi$ANNO == i,]
print(paste0("Bassi - Anno 20",i, sep=""))
print(gos_stat(x$DECESSI,i))
}
par(mfrow = c(3, 2))
for(i in 15:20) {
x <- df_aggr_medi[df_aggr_medi$ANNO == i,]
print(paste0("Medi - Anno 20",i, sep=""))
print(gos_stat(x$DECESSI,i))
}
par(mfrow = c(3, 2))
for(i in 15:20) {
x <- df_aggr_alti[df_aggr_alti$ANNO == i & !(df_aggr_alti$NOME_COMUNE %in% c("Milano","Bergamo","Brescia")),]
print(paste0("Alti - Anno 20",i, sep=""))
print(gos_stat(x$DECESSI,i))
}
$$
par(mfrow = c(3, 2))
for(i in 15:20) {
x <- df_aggr_medi[df_aggr_medi$ANNO == i,]
print(paste0("Medi - Anno 20",i, sep=""))
print(gos_stat(x$DECESSI,i))
}
## [1] "Medi - Anno 2015"
## Fitting of the distribution ' lnorm ' by maximum likelihood
## Parameters:
## estimate Std. Error
## meanlog 2.2537817 0.06361485
## sdlog 0.7225258 0.04498210
## Goodness-of-fit statistics
## Weibull Lognormal Gamma Loglogistic
## Kolmogorov-Smirnov statistic 0.1065330 0.09847301 0.09556002 0.06917855
## Cramer-von Mises statistic 0.2891741 0.10406118 0.18862769 0.08570355
## Anderson-Darling statistic 1.6267770 0.88629287 1.09581364 0.66380876
## Pareto
## Kolmogorov-Smirnov statistic 0.2319129
## Cramer-von Mises statistic 1.4088299
## Anderson-Darling statistic 7.8402495
##
## Goodness-of-fit criteria
## Weibull Lognormal Gamma Loglogistic Pareto
## Akaike's Information Criterion 867.2280 867.7113 862.7068 865.0353 903.6055
## Bayesian Information Criterion 872.9476 873.4309 868.4264 870.7549 909.3251
## [1] "Medi - Anno 2016"
## Fitting of the distribution ' lnorm ' by maximum likelihood
## Parameters:
## estimate Std. Error
## meanlog 2.1438531 0.06529743
## sdlog 0.7445053 0.04617188
## Goodness-of-fit statistics
## Weibull Lognormal Gamma Loglogistic
## Kolmogorov-Smirnov statistic 0.07767132 0.09011304 0.05849273 0.07211654
## Cramer-von Mises statistic 0.10900192 0.14809059 0.06266015 0.09205998
## Anderson-Darling statistic 0.77918414 0.96423074 0.48905673 0.78930128
## Pareto
## Kolmogorov-Smirnov statistic 0.1753318
## Cramer-von Mises statistic 1.1337252
## Anderson-Darling statistic 6.5056750
##
## Goodness-of-fit criteria
## Weibull Lognormal Gamma Loglogistic Pareto
## Akaike's Information Criterion 855.3264 853.6166 850.8753 856.0882 886.5424
## Bayesian Information Criterion 861.0615 859.3517 856.6104 861.8232 892.2775
## [1] "Medi - Anno 2017"
## Fitting of the distribution ' lnorm ' by maximum likelihood
## Parameters:
## estimate Std. Error
## meanlog 2.1692309 0.06917373
## sdlog 0.7887019 0.04891286
## Goodness-of-fit statistics
## Weibull Lognormal Gamma Loglogistic
## Kolmogorov-Smirnov statistic 0.1059538 0.07482570 0.08889435 0.06783799
## Cramer-von Mises statistic 0.1854201 0.08989504 0.12328535 0.08239378
## Anderson-Darling statistic 1.1421763 0.65550535 0.77660169 0.60159013
## Pareto
## Kolmogorov-Smirnov statistic 0.1822968
## Cramer-von Mises statistic 0.8740724
## Anderson-Darling statistic 5.2279817
##
## Goodness-of-fit criteria
## Weibull Lognormal Gamma Loglogistic Pareto
## Akaike's Information Criterion 876.6234 875.2087 872.6376 876.3367 900.5708
## Bayesian Information Criterion 882.3585 880.9437 878.3727 882.0718 906.3058
## [1] "Medi - Anno 2018"
## Fitting of the distribution ' lnorm ' by maximum likelihood
## Parameters:
## estimate Std. Error
## meanlog 2.2371787 0.06138730
## sdlog 0.6999229 0.04340698
## Goodness-of-fit statistics
## Weibull Lognormal Gamma Loglogistic
## Kolmogorov-Smirnov statistic 0.1205721 0.07376718 0.1131095 0.07981057
## Cramer-von Mises statistic 0.3592205 0.10822069 0.2780645 0.12439962
## Anderson-Darling statistic 2.0037181 0.63073832 1.4401341 0.73988129
## Pareto
## Kolmogorov-Smirnov statistic 0.2210531
## Cramer-von Mises statistic 1.2986285
## Anderson-Darling statistic 7.5409603
##
## Goodness-of-fit criteria
## Weibull Lognormal Gamma Loglogistic Pareto
## Akaike's Information Criterion 874.1105 861.8263 866.4812 864.5163 906.8896
## Bayesian Information Criterion 879.8455 867.5614 872.2163 870.2514 912.6247
## [1] "Medi - Anno 2019"
## Fitting of the distribution ' lnorm ' by maximum likelihood
## Parameters:
## estimate Std. Error
## meanlog 2.1756562 0.07311562
## sdlog 0.8304339 0.05170022
## Goodness-of-fit statistics
## Weibull Lognormal Gamma Loglogistic
## Kolmogorov-Smirnov statistic 0.07941785 0.07714341 0.08283226 0.09222728
## Cramer-von Mises statistic 0.12695871 0.12215267 0.10805279 0.12538874
## Anderson-Darling statistic 0.86175047 0.88274052 0.75445614 0.99702174
## Pareto
## Kolmogorov-Smirnov statistic 0.1602894
## Cramer-von Mises statistic 0.5496101
## Anderson-Darling statistic 3.7537903
##
## Goodness-of-fit criteria
## Weibull Lognormal Gamma Loglogistic Pareto
## Akaike's Information Criterion 881.0664 883.4672 879.3193 889.3455 901.6016
## Bayesian Information Criterion 886.7860 889.1868 885.0389 895.0651 907.3212
## [1] "Medi - Anno 2020"
## Fitting of the distribution ' lnorm ' by maximum likelihood
## Parameters:
## estimate Std. Error
## meanlog 3.6764143 0.02364755
## sdlog 0.2696236 0.01672031
## Goodness-of-fit statistics
## Weibull Lognormal Gamma Loglogistic
## Kolmogorov-Smirnov statistic 0.1418805 0.1113501 0.1260236 0.0934762
## Cramer-von Mises statistic 0.5587452 0.2332009 0.3202473 0.1869218
## Anderson-Darling statistic 3.5173490 1.6167444 2.1178753 1.5650559
## Pareto
## Kolmogorov-Smirnov statistic 0.4699311
## Cramer-von Mises statistic 6.6977298
## Anderson-Darling statistic 32.1513442
##
## Goodness-of-fit criteria
## Weibull Lognormal Gamma Loglogistic Pareto
## Akaike's Information Criterion 1012.928 988.0023 993.1281 996.5303 1229.578
## Bayesian Information Criterion 1018.663 993.7374 998.8632 1002.2653 1235.313
Violata quindi la seconda ipotesi per l’analisi della varianza passiamo a verificare la terza, quella di omeoschedasticità. Notiamo però che il p-value del test di Bartlett è sempre \(< 0.05\). Rifiutiamo quindi l’ipotesi nulla di omeoschedasticità.
bartlett.test(df_aggr_bassi$DECESSI, df_aggr_bassi$ANNO)
##
## Bartlett test of homogeneity of variances
##
## data: df_aggr_bassi$DECESSI and df_aggr_bassi$ANNO
## Bartlett's K-squared = 250.78, df = 5, p-value < 2.2e-16
bartlett.test(df_aggr_medi$DECESSI, df_aggr_medi$ANNO)
##
## Bartlett test of homogeneity of variances
##
## data: df_aggr_medi$DECESSI and df_aggr_medi$ANNO
## Bartlett's K-squared = 28.332, df = 5, p-value = 3.135e-05
bartlett.test(df_aggr_alti$DECESSI, df_aggr_alti$ANNO)
##
## Bartlett test of homogeneity of variances
##
## data: df_aggr_alti$DECESSI and df_aggr_alti$ANNO
## Bartlett's K-squared = 19.659, df = 5, p-value = 0.001448
Essendo le due ipotesi per l’ANOVA violate si può procedere con il test non parametrico di Friedman, e con il post-hoc test di Wilcoxon.
Il test ci mostra che c’è almeno una distribuzione differente dalle altre per ogni classe di mortalità, p-value motlo inferiore anche a \(0.01\).
## ANOVA per distribuzioni non normali ed eteroschedastiche, con varianze non omogenee,
## Procediamo con il test di Friedman, e con il post-hoc test di Wilcoxon
# https://www.r-statistics.com/2010/02/post-hoc-analysis-for-friedmans-test-r-code/
df_aggr_bassi <- df_aggr_bassi %>% rstatix::convert_as_factor(ANNO)
df_aggr_bassi %>% count(NOME_COMUNE) %>% filter(n < 6) %>% dplyr::select(NOME_COMUNE) -> comuni_del
df_aggr_bassi_Friedman <- df_aggr_bassi[!(df_aggr_bassi$NOME_COMUNE %in% comuni_del$NOME_COMUNE),]
df_aggr_bassi_15 <- df_aggr_bassi_Friedman[df_aggr_bassi_Friedman$ANNO == 15,]$DECESSI
df_aggr_bassi_16 <- df_aggr_bassi_Friedman[df_aggr_bassi_Friedman$ANNO == 16,]$DECESSI
df_aggr_bassi_17 <- df_aggr_bassi_Friedman[df_aggr_bassi_Friedman$ANNO == 17,]$DECESSI
df_aggr_bassi_18 <- df_aggr_bassi_Friedman[df_aggr_bassi_Friedman$ANNO == 18,]$DECESSI
df_aggr_bassi_19 <- df_aggr_bassi_Friedman[df_aggr_bassi_Friedman$ANNO == 19,]$DECESSI
df_aggr_bassi_20 <- df_aggr_bassi_Friedman[df_aggr_bassi_Friedman$ANNO == 20,]$DECESSI
df_aggr_bassi_matrix <- cbind(df_aggr_bassi_15, df_aggr_bassi_16, df_aggr_bassi_17,df_aggr_bassi_18,df_aggr_bassi_19,df_aggr_bassi_20)
#<0.05 differenza significativa tra almeno uno dei gruppi
friedman.test(df_aggr_bassi_matrix)
##
## Friedman rank sum test
##
## data: df_aggr_bassi_matrix
## Friedman chi-squared = 629.59, df = 5, p-value < 2.2e-16
#Medi
df_aggr_medi <- df_aggr_medi %>% rstatix::convert_as_factor(ANNO)
df_aggr_medi %>% count(NOME_COMUNE) %>% filter(n < 6) %>% dplyr::select(NOME_COMUNE) -> comuni_del
df_aggr_medi_Friedman <- df_aggr_medi[!(df_aggr_medi$NOME_COMUNE %in% comuni_del$NOME_COMUNE),]
df_aggr_medi_15 <- df_aggr_medi_Friedman[df_aggr_medi_Friedman$ANNO == 15,]$DECESSI
df_aggr_medi_16 <- df_aggr_medi_Friedman[df_aggr_medi_Friedman$ANNO == 16,]$DECESSI
df_aggr_medi_17 <- df_aggr_medi_Friedman[df_aggr_medi_Friedman$ANNO == 17,]$DECESSI
df_aggr_medi_18 <- df_aggr_medi_Friedman[df_aggr_medi_Friedman$ANNO == 18,]$DECESSI
df_aggr_medi_19 <- df_aggr_medi_Friedman[df_aggr_medi_Friedman$ANNO == 19,]$DECESSI
df_aggr_medi_20 <- df_aggr_medi_Friedman[df_aggr_medi_Friedman$ANNO == 20,]$DECESSI
df_aggr_medi_matrix <- cbind(df_aggr_medi_15, df_aggr_medi_16, df_aggr_medi_17,df_aggr_medi_18,df_aggr_medi_19,df_aggr_medi_20)
friedman.test(df_aggr_medi_matrix)
##
## Friedman rank sum test
##
## data: df_aggr_medi_matrix
## Friedman chi-squared = 282.46, df = 5, p-value < 2.2e-16
#Alti
df_aggr_alti <- df_aggr_alti %>% rstatix::convert_as_factor(ANNO)
df_aggr_alti %>% count(NOME_COMUNE) %>% filter(n < 6) %>% dplyr::select(NOME_COMUNE) -> comuni_del
df_aggr_alti_Friedman <- df_aggr_alti[!(df_aggr_alti$NOME_COMUNE %in% comuni_del$NOME_COMUNE),]
df_aggr_alti_15 <- df_aggr_alti_Friedman[df_aggr_alti_Friedman$ANNO == 15,]$DECESSI
df_aggr_alti_16 <- df_aggr_alti_Friedman[df_aggr_alti_Friedman$ANNO == 16,]$DECESSI
df_aggr_alti_17 <- df_aggr_alti_Friedman[df_aggr_alti_Friedman$ANNO == 17,]$DECESSI
df_aggr_alti_18 <- df_aggr_alti_Friedman[df_aggr_alti_Friedman$ANNO == 18,]$DECESSI
df_aggr_alti_19 <- df_aggr_alti_Friedman[df_aggr_alti_Friedman$ANNO == 19,]$DECESSI
df_aggr_alti_20 <- df_aggr_alti_Friedman[df_aggr_alti_Friedman$ANNO == 20,]$DECESSI
block <- 1:length(df_aggr_alti_20)
df_aggr_alti_matrix <- cbind(df_aggr_alti_15, df_aggr_alti_16, df_aggr_alti_17,df_aggr_alti_18,df_aggr_alti_19,df_aggr_alti_20,block)
friedman.test(df_aggr_alti_matrix)
##
## Friedman rank sum test
##
## data: df_aggr_alti_matrix
## Friedman chi-squared = 106.6, df = 6, p-value < 2.2e-16
Analizzando con il test post-hoc di Wilcox la differenza tra le distribuzioni a due a due notiamo che il \(2020\) è statisticamente differente da ogni altro anno.
Abbiamo dei valori significativi anche per altre distribuzioni, come per esempio nei comuni ad alta mortalità il \(2017\) è significativamente differente da tutti gli altri anni, come avevamo già notato precedentemente dall’analisi grafica.
#<0.05 differenza significativa
pairwise.wilcox.test(df_aggr_bassi_Friedman$DECESSI,df_aggr_bassi_Friedman$ANNO, paired = TRUE, exact=FALSE, p.adj="none", correct=FALSE)
##
## Pairwise comparisons using Wilcoxon signed rank test
##
## data: df_aggr_bassi_Friedman$DECESSI and df_aggr_bassi_Friedman$ANNO
##
## 15 16 17 18 19
## 16 0.020 - - - -
## 17 0.115 0.798 - - -
## 18 0.692 0.119 0.160 - -
## 19 0.782 0.013 0.043 0.685 -
## 20 <2e-16 <2e-16 <2e-16 <2e-16 <2e-16
##
## P value adjustment method: none
pairwise.wilcox.test(df_aggr_medi_Friedman$DECESSI,df_aggr_medi_Friedman$ANNO, paired = TRUE, exact=FALSE, p.adj="none", correct=FALSE)
##
## Pairwise comparisons using Wilcoxon signed rank test
##
## data: df_aggr_medi_Friedman$DECESSI and df_aggr_medi_Friedman$ANNO
##
## 15 16 17 18 19
## 16 0.048 - - - -
## 17 0.327 0.263 - - -
## 18 0.636 0.178 0.269 - -
## 19 0.434 0.101 0.701 0.743 -
## 20 <2e-16 <2e-16 <2e-16 <2e-16 <2e-16
##
## P value adjustment method: none
pairwise.wilcox.test(df_aggr_alti_Friedman$DECESSI,df_aggr_alti_Friedman$ANNO, paired = TRUE, exact=FALSE, p.adj="none", correct=FALSE)
##
## Pairwise comparisons using Wilcoxon signed rank test
##
## data: df_aggr_alti_Friedman$DECESSI and df_aggr_alti_Friedman$ANNO
##
## 15 16 17 18 19
## 16 0.6258 - - - -
## 17 0.0333 0.0357 - - -
## 18 0.4717 0.1580 0.0010 - -
## 19 0.1291 0.0591 0.0023 0.7826 -
## 20 5.2e-08 5.2e-08 5.2e-08 5.3e-08 5.2e-08
##
## P value adjustment method: none
Possiamo quindi concludere questa parte di analisi confermando che i comuni Lombardi, nel mese di Marzo \(2020\), hanno avuto un incremento significativo nel numero dei decessi.
# Prendo la Lombardia
df %>% filter(NOME_REGIONE == 'Lombardia')%>%
dplyr::select(NOME_PROVINCIA,NOME_COMUNE,DECESSI,CL_ETA,ANNO,MESE,GIORNO,DATA_INIZIO_DIFF, DATA, SESSO) -> df
#table(df$NOME_COMUNE) # alcuni comuni sono venuti con i nomi sbagliati li sistemo
df$NOME_COMUNE <- as.character(df$NOME_COMUNE)
df %>%
mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Alm\xe8", "Almè")) %>%
as.data.frame() -> df
# 2
df %>%
mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Barzan\xf2", "Barzanò"))%>%
as.data.frame() -> df
# 3
df %>%
mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Bascap\xe8", "Bascapè"))%>%
as.data.frame() -> df
# 7
df %>%
mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Cant\xf9", "Cantù"))%>%
as.data.frame() -> df
# 11
df %>%
mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Fenegr\xf2", "Fenegrò"))%>%
as.data.frame() -> df
df %>%
mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Gambol\xf2", "Gambolò"))%>%
as.data.frame() -> df
#
# 14
df %>%
mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Mont\xf9 Beccaria", "Montù Beccaria"))%>%
as.data.frame() -> df
# 15
df %>%
mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Muggi\xf2", "Muggiò"))%>%
as.data.frame() -> df
# 18
df %>%
mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Ro\xe8 Volciano", "Roè Volciano"))%>%
as.data.frame() -> df
# 19
df %>%
mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Sal\xf2", "Salò"))%>%
as.data.frame() -> df
df %>%
mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Santa Maria Ho\xe8", "Santa Maria Hoè"))%>%
as.data.frame() -> df
# 23
df %>%
mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Tem\xf9", "Temù"))%>%
as.data.frame() -> df
# 25
df %>%
mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Travac\xf2 Siccomario", "Travacò Siccomario"))%>%
as.data.frame() -> df
# 28
df %>%
mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Vigan\xf2", "Viganò"))%>%
as.data.frame() -> df
# 29
df %>%
mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Viggi\xf9", "Viggiù"))%>%
as.data.frame() -> df
# 30
df %>%
mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Villa d'Alm\xe8", "Villa d'Almè"))%>%
as.data.frame() -> df
# 31
df %>%
mutate(NOME_COMUNE=replace(NOME_COMUNE, NOME_COMUNE=="Zerbol\xf2", "Zerbolò"))%>%
as.data.frame() -> df
df_pop <- df_pop[-c(1,3:12,14,15)]
names(df_pop)[2] <- "POPOLAZIONE"
names(df_pop)[1] <- "NOME_COMUNE"
df %>% left_join(df_pop, by='NOME_COMUNE') -> df
## Warning: Column `NOME_COMUNE` joining character vector and factor, coercing into
## character vector
summary(df)
# controllo NA in popolazione
y<-dplyr::filter(df,is.na(POPOLAZIONE))
table(y$NOME_COMUNE)
Decido di sistemare Manualmente per questi comuni la popolazione totale
df %>%
mutate(POPOLAZIONE=replace(POPOLAZIONE, NOME_COMUNE=="Borgo Mantovano", 5529))%>%
as.data.frame() -> df
df %>%
mutate(POPOLAZIONE=replace(POPOLAZIONE, NOME_COMUNE=="Borgocarbonara", 1931))%>%
as.data.frame() -> df
df %>%
mutate(POPOLAZIONE=replace(POPOLAZIONE, NOME_COMUNE=="Cadrezzate con Osmate", 2659))%>%
as.data.frame() -> df
df %>%
mutate(POPOLAZIONE=replace(POPOLAZIONE, NOME_COMUNE=="Castelgerundo", 1498))%>%
as.data.frame() -> df
df %>%
mutate(POPOLAZIONE=replace(POPOLAZIONE, NOME_COMUNE=="Centro Valle Intelvi", 3227))%>%
as.data.frame() -> df
df %>%
mutate(POPOLAZIONE=replace(POPOLAZIONE, NOME_COMUNE=="Colli Verdi", 1074))%>%
as.data.frame() -> df
df %>%
mutate(POPOLAZIONE=replace(POPOLAZIONE, NOME_COMUNE=="Gornate Olona", 2223))%>%
as.data.frame() -> df
df %>%
mutate(POPOLAZIONE=replace(POPOLAZIONE, NOME_COMUNE=="Piadena Drizzona", 563))%>%
as.data.frame() -> df
df %>%
mutate(POPOLAZIONE=replace(POPOLAZIONE, NOME_COMUNE=="Puegnago del Garda", 3439))%>%
as.data.frame() -> df
df %>%
mutate(POPOLAZIONE=replace(POPOLAZIONE, NOME_COMUNE=="San Giorgio Bigarello", 13654))%>%
as.data.frame() -> df
df %>%
mutate(POPOLAZIONE=replace(POPOLAZIONE, NOME_COMUNE=="Solbiate con Cagno", 4653))%>%
as.data.frame() -> df
df %>%
mutate(POPOLAZIONE=replace(POPOLAZIONE, NOME_COMUNE=="Torre de' Busi", 2130))%>%
as.data.frame() -> df
df %>%
mutate(POPOLAZIONE=replace(POPOLAZIONE, NOME_COMUNE=="Valvarrone", 551))%>%
as.data.frame() -> df
df %>%
mutate(POPOLAZIONE=replace(POPOLAZIONE, NOME_COMUNE=="Vermezzo con Zelo", 5805))%>%
as.data.frame() -> df
Sono presenti però anche i dati per le provincie(max 3250315 che equivale alla provincia di milano,) decido di sistemare
# Brescia
df %>% filter(POPOLAZIONE!=1265954) ->df
# milano
df %>% filter(POPOLAZIONE!=3250315) -> df
# bergamo
df %>% filter(POPOLAZIONE!=1114590) -> df
# como
df %>% filter(POPOLAZIONE!=599204) -> df
# cremona
df %>% filter(POPOLAZIONE!=358955) -> df
# lecco
df %>% filter(POPOLAZIONE!=337380) -> df
#lodi
df %>% filter(POPOLAZIONE!=230198) -> df
# mantova
df %>% filter(POPOLAZIONE!=412292) -> df
#Pavia
df %>% filter(POPOLAZIONE!=545888) -> df
# sondrio
df %>% filter(POPOLAZIONE!=181095) -> df
# varese
df %>% filter(POPOLAZIONE!=890768) -> df
summary(df)
## NOME_PROVINCIA NOME_COMUNE DECESSI
## Milano :344424 Length:1762884 Min. : 0.0
## Brescia :240672 Class :character 1st Qu.: 0.0
## Bergamo :238272 Mode :character Median : 0.0
## Varese :173100 Mean : 712.7
## Monza e della Brianza:141540 3rd Qu.: 0.0
## Pavia :140304 Max. :9999.0
## (Other) :484572
## CL_ETA ANNO MESE GIORNO
## Min. : 0.00 Length:1762884 Min. :1.000 Min. : 1.00
## 1st Qu.:15.00 Class :character 1st Qu.:1.000 1st Qu.: 8.00
## Median :17.00 Mode :character Median :2.000 Median :16.00
## Mean :16.54 Mean :2.415 Mean :15.55
## 3rd Qu.:18.00 3rd Qu.:3.000 3rd Qu.:23.00
## Max. :21.00 Max. :4.000 Max. :31.00
##
## DATA_INIZIO_DIFF DATA SESSO
## 1 aprile :906396 Min. :2015-01-01 Length:1762884
## 16 aprile : 99204 1st Qu.:2016-02-28 Class :character
## 8 aprile :195672 Median :2017-08-31 Mode :character
## Dati 2020 n.d.:561612 Mean :2017-08-28
## 3rd Qu.:2019-02-28
## Max. :2020-04-30
## NA's :3400
## POPOLAZIONE
## Min. : 33
## 1st Qu.: 4355
## Median : 8883
## Mean : 34510
## 3rd Qu.: 20625
## Max. :1378689
##
y<-dplyr::filter(df,is.na(DATA))
df %>% anti_join(y, by="DATA") -> df
summary(df)
## NOME_PROVINCIA NOME_COMUNE DECESSI
## Milano :343728 Length:1759484 Min. : 0.0
## Brescia :240200 Class :character 1st Qu.: 0.0
## Bergamo :237768 Mode :character Median : 0.0
## Varese :172812 Mean : 714.1
## Monza e della Brianza:141260 3rd Qu.: 0.0
## Pavia :140096 Max. :9999.0
## (Other) :483620
## CL_ETA ANNO MESE GIORNO
## Min. : 0.00 Length:1759484 Min. :1.000 Min. : 1.00
## 1st Qu.:15.00 Class :character 1st Qu.:1.000 1st Qu.: 8.00
## Median :17.00 Mode :character Median :2.000 Median :16.00
## Mean :16.54 Mean :2.416 Mean :15.52
## 3rd Qu.:18.00 3rd Qu.:3.000 3rd Qu.:23.00
## Max. :21.00 Max. :4.000 Max. :31.00
##
## DATA_INIZIO_DIFF DATA SESSO
## 1 aprile :904396 Min. :2015-01-01 Length:1759484
## 16 aprile : 99060 1st Qu.:2016-02-28 Class :character
## 8 aprile :195192 Median :2017-08-31 Mode :character
## Dati 2020 n.d.:560836 Mean :2017-08-28
## 3rd Qu.:2019-02-28
## Max. :2020-04-30
##
## POPOLAZIONE
## Min. : 33
## 1st Qu.: 4355
## Median : 8883
## Mean : 34469
## 3rd Qu.: 20544
## Max. :1378689
##
df <- subset(df, !(GIORNO==29 & MESE==2))
df <- subset(df , MESE!=4)
Per comodità eliminiamo tutti
# metto a NA i valori 9999
df[df == 9999] <- NA
df %>% filter(POPOLAZIONE>=60000) -> df_60000
y<-dplyr::filter(df_60000,is.na(DECESSI))
table(y$NOME_COMUNE)
##
## Vigevano
## 1116
# DEVO TOGLIERE VIGEVANO POICHÈ HO 1116 OSSERVAZIONI A NA
df_60000 %>% filter(NOME_COMUNE!='Vigevano') -> df_60000
table(df_60000$NOME_COMUNE)
##
## Bergamo Brescia Busto Arsizio Cinisello Balsamo
## 9240 10740 7524 7260
## Como Cremona Legnano Milano
## 8268 8064 6348 15552
## Monza Pavia Sesto San Giovanni Varese
## 8796 7764 7788 8040
Abbiamno deciso di mostrare le serie storiche di Comuni con più di 60.000 abitanti:
Inizialmente prendiamo il dato assoluto dei decessi per cui sostanzialmente andiamo a fare un ranking delle morti con un incremento percentuale maggiore.
listacomuni <- c("Bergamo", "Brescia", "Busto Arsizio", "Cinisello Balsamo", "Como", "Cremona", "Legnano", "Milano", "Monza", "Pavia", "Sesto San Giovanni", "Varese")
datasets <- list()
for(i in listacomuni){
df_60000 %>% filter(NOME_COMUNE==i) %>%
mutate(CL_ETA= case_when(CL_ETA<= 6 ~ "GIOVANI",
CL_ETA %in% c(7,8,9,10,11,12,13,14) ~ "MEDI",
CL_ETA > 14 ~ "ANZIANI",)) %>%
mutate(DECESSI=replace(DECESSI, DECESSI==9999, NA)) -> x
x %>% group_by(DATA) %>%
summarise(DECESSI = sum(DECESSI)) %>%
as_tibble() -> x
assign(paste0('df_comune_',i), x)
datasets <- append(datasets,paste0('df_comune_',i))
}
# Milano
ts_Bergamo <- zoo(df_comune_Bergamo$DECESSI, df_comune_Bergamo$DATA)
ts_Bergamo <- ts(ts_Bergamo, frequency = 90)
autoplot(ts_Bergamo, main = paste0('Decessi Bergamo 2015-2020'), xlab = 'Anni a 90 giorni', ylab = 'Decessi')
Il codice per il plot si ripete.
ts_total <- ts.union(ts_Bergamo, ts_Brescia, ts_Bs, ts_Cb, ts_Como, ts_Cremona, ts_Legnano, ts_milano, ts_Monza, ts_Pavia, ts_ssg, ts_Varese)
colnames(ts_total)
## [1] "ts_Bergamo" "ts_Brescia" "ts_Bs" "ts_Cb" "ts_Como"
## [6] "ts_Cremona" "ts_Legnano" "ts_milano" "ts_Monza" "ts_Pavia"
## [11] "ts_ssg" "ts_Varese"
colnames(ts_total)[1:12] <- listacomuni
autoplot(ts_total, main = paste0('Decessi Milano'), xlab = 'Marzo', ylab = 'Decessi')
Come possiamo notare Milano in termini assoluti ha avuto l’incremento maggiore per cui nel mese di Marzo ha avuto un rapido incremento fino ad arrivare per vari giorni a 80 decessi al giorno. Da notare anche Bergamo con la media di 30 decessi al giorno nel mese di Marzo.
Decidiamo quindi di creare una variabile che chiameremo “Decessi pro capite” per cui andiamo a stimare le morti giornaliere rispetto alla popolazione (stima del 2019) andando a vedere non più i termini assoluti. Vogliamo vedere se le zone più colpite rimangono le stesse.
listacomuni <- c("Bergamo", "Brescia", "Busto Arsizio", "Cinisello Balsamo", "Como", "Cremona", "Legnano", "Milano", "Monza", "Pavia", "Sesto San Giovanni", "Varese")
datasets <- list()
str(x$DECESSISUPOPOLAZIONE)
## Warning: Unknown or uninitialised column: `DECESSISUPOPOLAZIONE`.
## NULL
for(i in listacomuni){
df_60000 %>% filter(NOME_COMUNE==i) %>%
mutate(CL_ETA= case_when(CL_ETA<= 6 ~ "GIOVANI",
CL_ETA %in% c(7,8,9,10,11,12,13,14) ~ "MEDI",
CL_ETA > 14 ~ "ANZIANI",)) %>%
mutate(DECESSISUPOPOLAZIONE = (DECESSI/POPOLAZIONE)) %>%
mutate(DECESSI=replace(DECESSI, DECESSI==9999, NA)) -> x
x %>% group_by(DATA) %>%
summarise(DECESSISUPOPOLAZIONE = sum(DECESSISUPOPOLAZIONE)) %>%
as_tibble() -> x
assign(paste0('df_comune_',i), x)
datasets <- append(datasets,paste0('df_comune_',i))
}
# Milano
ts_Bergamo <- zoo(df_comune_Bergamo$DECESSISUPOPOLAZIONE, df_comune_Bergamo$DATA)
ts_Bergamo <- ts(ts_Bergamo, frequency = 90)
autoplot(ts_Bergamo, main = paste0('Decessi Bergamo 2015-2020'), xlab = 'Anni a 90 giorni', ylab = 'Decessi')
Il codice per il plot si ripete.
ts_total <- ts.union(ts_Bergamo, ts_Brescia, ts_Bs, ts_Cb, ts_Como, ts_Cremona, ts_Legnano, ts_milano, ts_Monza, ts_Pavia, ts_ssg, ts_Varese)
# colnames(ts_total)
colnames(ts_total)[1:12] <- listacomuni
autoplot(ts_total, main = paste0('Decessi Milano'), xlab = 'Marzo', ylab = 'Decessi')
Al contrario del grafico precedente se consideriamo i morti rispetto alla popolazione notiamo che non è Milano la città con il record più negativo, ma diventa Bergamo che è seguito dal Cremona e Brescia.
Decidiamo di analizzare i NULL, per cui secondo le specifiche del dataset questi comuni non hanno avuto un incremento del 20% della mortalità dell’anno 2020. Perciò decidiamo arbitrariamente di sostituirli con i valori delle medie dei 5 anni precedenti. Così facendo notiamo fin da subito un incremento molto evidente della mortalità per il mese di Marzo anno 2020. Supponiamo perciò si tratti non solo unicamente del Covid poichè è per noio impossibile da dimostrare, ma sicuramente c’è stato un forte impatto.
df %>% filter(ANNO %in% c(15,16,17,18,19),
MESE %in% c(1,2,3),
DECESSI != 9999) %>%
group_by(NOME_COMUNE,ANNO,MESE) %>%
summarise(DECESSI = sum(DECESSI)) %>% group_by(NOME_COMUNE,MESE) %>%
summarise(DECESSI=mean(DECESSI)) -> mean_dataset
df %>% filter(DECESSI == 9999,
MESE %in% c(1,2,3)) %>%
mutate(DECESSI = NA) %>%
group_by(NOME_COMUNE,ANNO,MESE) %>%
summarise(DECESSI = sum()) %>%
left_join(mean_dataset, by = c("NOME_COMUNE" = "NOME_COMUNE", "MESE" = "MESE")) %>% mutate(DECESSI.x = as.integer(DECESSI.y)) -> finale
df %>%filter(MESE %in% c(1,2,3)) %>% mutate(DECESSI = ifelse(DECESSI==9999, NA, DECESSI)) %>% group_by(NOME_COMUNE,ANNO,MESE) %>%
summarise(DECESSI = sum(DECESSI)) -> oprpr
finale_sistemati_null <- oprpr %>% left_join(finale, by = c("NOME_COMUNE" = "NOME_COMUNE", "MESE" = "MESE")) %>% mutate(DECESSI = ifelse(is.na(DECESSI), as.integer(DECESSI.y), DECESSI))
colnames(finale_sistemati_null)[2] <- "ANNO"
finale_sistemati_null <- finale_sistemati_null[-c(5:7)]
finale_sistemati_null %>%
group_by(ANNO, MESE) %>%
summarise(DECESSI = sum(DECESSI)) -> finale_sistemati_null